/*
 *    This program is free software; you can redistribute it and/or modify
 *    it under the terms of the GNU General Public License as published by
 *    the Free Software Foundation; either version 2 of the License, or (at
 *    your option) any later version.
 *
 *    This program is distributed in the hope that it will be useful, but
 *    WITHOUT ANY WARRANTY; without even the implied warranty of
 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 *    General Public License for more details.
 *
 *    You should have received a copy of the GNU General Public License
 *    along with this program; if not, write to the Free Software
 *    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.  */

/*
 *    PaceMatrix.java
 *    Copyright (C) 2002 University of Waikato, Hamilton, New Zealand
 *
 */

package weka.classifiers.functions.pace;

import weka.core.RevisionUtils;
import weka.core.matrix.DoubleVector;
import weka.core.matrix.FlexibleDecimalFormat;
import weka.core.matrix.IntVector;
import weka.core.matrix.Matrix;
import weka.core.matrix.Maths;

import java.util.Random;
import java.text.DecimalFormat;

/**
 * Class for matrix manipulation used for pace regression.
 * <p>
 * 
 * REFERENCES
 * <p>
 * 
 * Wang, Y. (2000). "A new approach to fitting linear models in high dimensional
 * spaces." PhD Thesis. Department of Computer Science, University of Waikato,
 * New Zealand.
 * <p>
 * 
 * Wang, Y. and Witten, I. H. (2002). "Modeling for optimal probability
 * prediction." Proceedings of ICML'2002. Sydney.
 * <p>
 * 
 * @author Yong Wang (yongwang@cs.waikato.ac.nz)
 * @version $Revision: 1.6 $
 */
public class PaceMatrix extends Matrix {

	/** for serialization */
	static final long serialVersionUID = 2699925616857843973L;

	/*
	 * ------------------------ Constructors ------------------------
	 */

	/**
	 * Construct an m-by-n PACE matrix of zeros.
	 * 
	 * @param m
	 *            Number of rows.
	 * @param n
	 *            Number of colums.
	 */
	public PaceMatrix(int m, int n) {
		super(m, n);
	}

	/**
	 * Construct an m-by-n constant PACE matrix.
	 * 
	 * @param m
	 *            Number of rows.
	 * @param n
	 *            Number of colums.
	 * @param s
	 *            Fill the matrix with this scalar value.
	 */
	public PaceMatrix(int m, int n, double s) {
		super(m, n, s);
	}

	/**
	 * Construct a PACE matrix from a 2-D array.
	 * 
	 * @param A
	 *            Two-dimensional array of doubles.
	 * @throws IllegalArgumentException
	 *             All rows must have the same length
	 */
	public PaceMatrix(double[][] A) {
		super(A);
	}

	/**
	 * Construct a PACE matrix quickly without checking arguments.
	 * 
	 * @param A
	 *            Two-dimensional array of doubles.
	 * @param m
	 *            Number of rows.
	 * @param n
	 *            Number of colums.
	 */
	public PaceMatrix(double[][] A, int m, int n) {
		super(A, m, n);
	}

	/**
	 * Construct a PaceMatrix from a one-dimensional packed array
	 * 
	 * @param vals
	 *            One-dimensional array of doubles, packed by columns (ala
	 *            Fortran).
	 * @param m
	 *            Number of rows.
	 * @throws IllegalArgumentException
	 *             Array length must be a multiple of m.
	 */
	public PaceMatrix(double vals[], int m) {
		super(vals, m);
	}

	/**
	 * Construct a PaceMatrix with a single column from a DoubleVector
	 * 
	 * @param v
	 *            DoubleVector
	 */
	public PaceMatrix(DoubleVector v) {
		this(v.size(), 1);
		setMatrix(0, v.size() - 1, 0, v);
	}

	/**
	 * Construct a PaceMatrix from a Matrix
	 * 
	 * @param X
	 *            Matrix
	 */
	public PaceMatrix(Matrix X) {
		super(X.getRowDimension(), X.getColumnDimension());
		A = X.getArray();
	}

	/*
	 * ------------------------ Public Methods ------------------------
	 */

	/**
	 * Set the row dimenion of the matrix
	 * 
	 * @param rowDimension
	 *            the row dimension
	 */
	public void setRowDimension(int rowDimension) {
		m = rowDimension;
	}

	/**
	 * Set the column dimenion of the matrix
	 * 
	 * @param columnDimension
	 *            the column dimension
	 */
	public void setColumnDimension(int columnDimension) {
		n = columnDimension;
	}

	/**
	 * Clone the PaceMatrix object.
	 * 
	 * @return the clone
	 */
	public Object clone() {
		PaceMatrix X = new PaceMatrix(m, n);
		double[][] C = X.getArray();
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				C[i][j] = A[i][j];
			}
		}
		return (Object) X;
	}

	/**
	 * Add a value to an element and reset the element
	 * 
	 * @param i
	 *            the row number of the element
	 * @param j
	 *            the column number of the element
	 * @param s
	 *            the double value to be added with
	 */
	public void setPlus(int i, int j, double s) {
		A[i][j] += s;
	}

	/**
	 * Multiply a value with an element and reset the element
	 * 
	 * @param i
	 *            the row number of the element
	 * @param j
	 *            the column number of the element
	 * @param s
	 *            the double value to be multiplied with
	 */
	public void setTimes(int i, int j, double s) {
		A[i][j] *= s;
	}

	/**
	 * Set the submatrix A[i0:i1][j0:j1] with a same value
	 * 
	 * @param i0
	 *            the index of the first element of the column
	 * @param i1
	 *            the index of the last element of the column
	 * @param j0
	 *            the index of the first column
	 * @param j1
	 *            the index of the last column
	 * @param s
	 *            the value to be set to
	 */
	public void setMatrix(int i0, int i1, int j0, int j1, double s) {
		try {
			for (int i = i0; i <= i1; i++) {
				for (int j = j0; j <= j1; j++) {
					A[i][j] = s;
				}
			}
		} catch (ArrayIndexOutOfBoundsException e) {
			throw new ArrayIndexOutOfBoundsException("Index out of bounds");
		}
	}

	/**
	 * Set the submatrix A[i0:i1][j] with the values stored in a DoubleVector
	 * 
	 * @param i0
	 *            the index of the first element of the column
	 * @param i1
	 *            the index of the last element of the column
	 * @param j
	 *            the index of the column
	 * @param v
	 *            the vector that stores the values
	 */
	public void setMatrix(int i0, int i1, int j, DoubleVector v) {
		for (int i = i0; i <= i1; i++) {
			A[i][j] = v.get(i - i0);
		}
	}

	/**
	 * Set the whole matrix from a 1-D array
	 * 
	 * @param v
	 *            1-D array of doubles
	 * @param columnFirst
	 *            Whether to fill the column first or the row.
	 * @throws ArrayIndexOutOfBoundsException
	 *             Submatrix indices
	 */
	public void setMatrix(double[] v, boolean columnFirst) {
		try {
			if (v.length != m * n)
				throw new IllegalArgumentException("sizes not match.");
			int i, j, count = 0;
			if (columnFirst) {
				for (i = 0; i < m; i++) {
					for (j = 0; j < n; j++) {
						A[i][j] = v[count];
						count++;
					}
				}
			} else {
				for (j = 0; j < n; j++) {
					for (i = 0; i < m; i++) {
						A[i][j] = v[count];
						count++;
					}
				}
			}

		} catch (ArrayIndexOutOfBoundsException e) {
			throw new ArrayIndexOutOfBoundsException("Submatrix indices");
		}
	}

	/**
	 * Returns the maximum absolute value of all elements
	 * 
	 * @return the maximum value
	 */
	public double maxAbs() {
		double ma = Math.abs(A[0][0]);
		for (int j = 0; j < n; j++) {
			for (int i = 0; i < m; i++) {
				ma = Math.max(ma, Math.abs(A[i][j]));
			}
		}
		return ma;
	}

	/**
	 * Returns the maximum absolute value of some elements of a column, that is,
	 * the elements of A[i0:i1][j].
	 * 
	 * @param i0
	 *            the index of the first element of the column
	 * @param i1
	 *            the index of the last element of the column
	 * @param j
	 *            the index of the column
	 * @return the maximum value
	 */
	public double maxAbs(int i0, int i1, int j) {
		double m = Math.abs(A[i0][j]);
		for (int i = i0 + 1; i <= i1; i++) {
			m = Math.max(m, Math.abs(A[i][j]));
		}
		return m;
	}

	/**
	 * Returns the minimum absolute value of some elements of a column, that is,
	 * the elements of A[i0:i1][j].
	 * 
	 * @param i0
	 *            the index of the first element of the column
	 * @param i1
	 *            the index of the last element of the column
	 * @param column
	 *            the index of the column
	 * @return the minimum value
	 */
	public double minAbs(int i0, int i1, int column) {
		double m = Math.abs(A[i0][column]);
		for (int i = i0 + 1; i <= i1; i++) {
			m = Math.min(m, Math.abs(A[i][column]));
		}
		return m;
	}

	/**
	 * Check if the matrix is empty
	 * 
	 * @return true if the matrix is empty
	 */
	public boolean isEmpty() {
		if (m == 0 || n == 0)
			return true;
		if (A == null)
			return true;
		return false;
	}

	/**
	 * Return a DoubleVector that stores a column of the matrix
	 * 
	 * @param j
	 *            the index of the column
	 * @return the column
	 */
	public DoubleVector getColumn(int j) {
		DoubleVector v = new DoubleVector(m);
		double[] a = v.getArray();
		for (int i = 0; i < m; i++)
			a[i] = A[i][j];
		return v;
	}

	/**
	 * Return a DoubleVector that stores some elements of a column of the matrix
	 * 
	 * @param i0
	 *            the index of the first element of the column
	 * @param i1
	 *            the index of the last element of the column
	 * @param j
	 *            the index of the column
	 * @return the DoubleVector
	 */
	public DoubleVector getColumn(int i0, int i1, int j) {
		DoubleVector v = new DoubleVector(i1 - i0 + 1);
		double[] a = v.getArray();
		int count = 0;
		for (int i = i0; i <= i1; i++) {
			a[count] = A[i][j];
			count++;
		}
		return v;
	}

	/**
	 * Multiplication between a row (or part of a row) of the first matrix and a
	 * column (or part or a column) of the second matrix.
	 * 
	 * @param i
	 *            the index of the row in the first matrix
	 * @param j0
	 *            the index of the first column in the first matrix
	 * @param j1
	 *            the index of the last column in the first matrix
	 * @param B
	 *            the second matrix
	 * @param l
	 *            the index of the column in the second matrix
	 * @return the result of the multiplication
	 */
	public double times(int i, int j0, int j1, PaceMatrix B, int l) {
		double s = 0.0;
		for (int j = j0; j <= j1; j++) {
			s += A[i][j] * B.A[j][l];
		}
		return s;
	}

	/**
	 * Decimal format for converting a matrix into a string
	 * 
	 * @return the default decimal format
	 */
	protected DecimalFormat[] format() {
		return format(0, m - 1, 0, n - 1, 7, false);
	}

	/**
	 * Decimal format for converting a matrix into a string
	 * 
	 * @param digits
	 *            the number of digits
	 * @return the decimal format
	 */
	protected DecimalFormat[] format(int digits) {
		return format(0, m - 1, 0, n - 1, digits, false);
	}

	/**
	 * Decimal format for converting a matrix into a string
	 * 
	 * @param digits
	 *            the number of digits
	 * @param trailing
	 * @return the decimal format
	 */
	protected DecimalFormat[] format(int digits, boolean trailing) {
		return format(0, m - 1, 0, n - 1, digits, trailing);
	}

	/**
	 * Decimal format for converting a matrix into a string
	 * 
	 * @param i0
	 * @param i1
	 * @param j
	 * @param digits
	 *            the number of digits
	 * @param trailing
	 * @return the decimal format
	 */
	protected DecimalFormat format(int i0, int i1, int j, int digits,
			boolean trailing) {
		FlexibleDecimalFormat df = new FlexibleDecimalFormat(digits, trailing);
		df.grouping(true);
		for (int i = i0; i <= i1; i++)
			df.update(A[i][j]);
		return df;
	}

	/**
	 * Decimal format for converting a matrix into a string
	 * 
	 * @param i0
	 * @param i1
	 * @param j0
	 * @param j1
	 * @param trailing
	 * @param digits
	 *            the number of digits
	 * @return the decimal format
	 */
	protected DecimalFormat[] format(int i0, int i1, int j0, int j1,
			int digits, boolean trailing) {
		DecimalFormat[] f = new DecimalFormat[j1 - j0 + 1];
		for (int j = j0; j <= j1; j++) {
			f[j] = format(i0, i1, j, digits, trailing);
		}
		return f;
	}

	/**
	 * Converts matrix to string
	 * 
	 * @return the matrix as string
	 */
	public String toString() {
		return toString(5, false);
	}

	/**
	 * Converts matrix to string
	 * 
	 * @param digits
	 *            number of digits after decimal point
	 * @param trailing
	 *            true if trailing zeros are padded
	 * @return the matrix as string
	 */
	public String toString(int digits, boolean trailing) {

		if (isEmpty())
			return "null matrix";

		StringBuffer text = new StringBuffer();
		DecimalFormat[] nf = format(digits, trailing);
		int numCols = 0;
		int count = 0;
		int width = 80;
		int lenNumber;

		int[] nCols = new int[n];
		int nk = 0;
		for (int j = 0; j < n; j++) {
			lenNumber = nf[j].format(A[0][j]).length();
			if (count + 1 + lenNumber > width - 1) {
				nCols[nk++] = numCols;
				count = 0;
				numCols = 0;
			}
			count += 1 + lenNumber;
			++numCols;
		}
		nCols[nk] = numCols;

		nk = 0;
		for (int k = 0; k < n;) {
			for (int i = 0; i < m; i++) {
				for (int j = k; j < k + nCols[nk]; j++)
					text.append(" " + nf[j].format(A[i][j]));
				text.append("\n");
			}
			k += nCols[nk];
			++nk;
			text.append("\n");
		}

		return text.toString();
	}

	/**
	 * Squared sum of a column or row in a matrix
	 * 
	 * @param j
	 *            the index of the column or row
	 * @param i0
	 *            the index of the first element
	 * @param i1
	 *            the index of the last element
	 * @param col
	 *            if true, sum over a column; otherwise, over a row
	 * @return the squared sum
	 */
	public double sum2(int j, int i0, int i1, boolean col) {
		double s2 = 0;
		if (col) { // column
			for (int i = i0; i <= i1; i++)
				s2 += A[i][j] * A[i][j];
		} else {
			for (int i = i0; i <= i1; i++)
				s2 += A[j][i] * A[j][i];
		}
		return s2;
	}

	/**
	 * Squared sum of columns or rows of a matrix
	 * 
	 * @param col
	 *            if true, sum over columns; otherwise, over rows
	 * @return the squared sum
	 */
	public double[] sum2(boolean col) {
		int l = col ? n : m;
		int p = col ? m : n;
		double[] s2 = new double[l];
		for (int i = 0; i < l; i++)
			s2[i] = sum2(i, 0, p - 1, col);
		return s2;
	}

	/**
	 * Constructs single Householder transformation for a column
	 * 
	 * @param j
	 *            the index of the column
	 * @param k
	 *            the index of the row
	 * @return d and q
	 */
	public double[] h1(int j, int k) {
		double dq[] = new double[2];
		double s2 = sum2(j, k, m - 1, true);
		dq[0] = A[k][j] >= 0 ? -Math.sqrt(s2) : Math.sqrt(s2);
		A[k][j] -= dq[0];
		dq[1] = A[k][j] * dq[0];
		return dq;
	}

	/**
	 * Performs single Householder transformation on one column of a matrix
	 * 
	 * @param j
	 *            the index of the column
	 * @param k
	 *            the index of the row
	 * @param q
	 *            q = - u'u/2; must be negative
	 * @param b
	 *            the matrix to be transformed
	 * @param l
	 *            the column of the matrix b
	 */
	public void h2(int j, int k, double q, PaceMatrix b, int l) {
		double s = 0, alpha;
		for (int i = k; i < m; i++)
			s += A[i][j] * b.A[i][l];
		alpha = s / q;
		for (int i = k; i < m; i++)
			b.A[i][l] += alpha * A[i][j];
	}

	/**
	 * Constructs the Givens rotation
	 * 
	 * @param a
	 * @param b
	 * @return a double array that stores the cosine and sine values
	 */
	public double[] g1(double a, double b) {
		double cs[] = new double[2];
		double r = Maths.hypot(a, b);
		if (r == 0.0) {
			cs[0] = 1;
			cs[1] = 0;
		} else {
			cs[0] = a / r;
			cs[1] = b / r;
		}
		return cs;
	}

	/**
	 * Performs the Givens rotation
	 * 
	 * @param cs
	 *            a array storing the cosine and sine values
	 * @param i0
	 *            the index of the row of the first element
	 * @param i1
	 *            the index of the row of the second element
	 * @param j
	 *            the index of the column
	 */
	public void g2(double cs[], int i0, int i1, int j) {
		double w = cs[0] * A[i0][j] + cs[1] * A[i1][j];
		A[i1][j] = -cs[1] * A[i0][j] + cs[0] * A[i1][j];
		A[i0][j] = w;
	}

	/**
	 * Forward ordering of columns in terms of response explanation. On input,
	 * matrices A and b are already QR-transformed. The indices of transformed
	 * columns are stored in the pivoting vector.
	 * 
	 * @param b
	 *            the PaceMatrix b
	 * @param pvt
	 *            the pivoting vector
	 * @param k0
	 *            the first k0 columns (in pvt) of A are not to be changed
	 **/
	public void forward(PaceMatrix b, IntVector pvt, int k0) {
		for (int j = k0; j < Math.min(pvt.size(), m); j++) {
			steplsqr(b, pvt, j, mostExplainingColumn(b, pvt, j), true);
		}
	}

	/**
	 * Returns the index of the column that has the largest (squared) response,
	 * when each of columns pvt[ks:] is moved to become the ks-th column. On
	 * input, A and b are both QR-transformed.
	 * 
	 * @param b
	 *            response
	 * @param pvt
	 *            pivoting index of A
	 * @param ks
	 *            columns pvt[ks:] of A are to be tested
	 * @return the index of the column
	 */
	public int mostExplainingColumn(PaceMatrix b, IntVector pvt, int ks) {
		double val;
		int[] p = pvt.getArray();
		double ma = columnResponseExplanation(b, pvt, ks, ks);
		int jma = ks;
		for (int i = ks + 1; i < pvt.size(); i++) {
			val = columnResponseExplanation(b, pvt, i, ks);
			if (val > ma) {
				ma = val;
				jma = i;
			}
		}
		return jma;
	}

	/**
	 * Backward ordering of columns in terms of response explanation. On input,
	 * matrices A and b are already QR-transformed. The indices of transformed
	 * columns are stored in the pivoting vector.
	 * 
	 * A and b must have the same number of rows, being the (pseudo-)rank.
	 * 
	 * @param b
	 *            PaceMatrix b
	 * @param pvt
	 *            pivoting vector
	 * @param ks
	 *            number of QR-transformed columns; psuedo-rank of A
	 * @param k0
	 *            first k0 columns in pvt[] are not to be ordered.
	 */
	public void backward(PaceMatrix b, IntVector pvt, int ks, int k0) {
		for (int j = ks; j > k0; j--) {
			steplsqr(b, pvt, j, leastExplainingColumn(b, pvt, j, k0), false);
		}
	}

	/**
	 * Returns the index of the column that has the smallest (squared) response,
	 * when the column is moved to become the (ks-1)-th column. On input, A and
	 * b are both QR-transformed.
	 * 
	 * @param b
	 *            response
	 * @param pvt
	 *            pivoting index of A
	 * @param ks
	 *            psudo-rank of A
	 * @param k0
	 *            A[][pvt[0:(k0-1)]] are excluded from the testing.
	 * @return the index of the column
	 */
	public int leastExplainingColumn(PaceMatrix b, IntVector pvt, int ks, int k0) {
		double val;
		int[] p = pvt.getArray();
		double mi = columnResponseExplanation(b, pvt, ks - 1, ks);
		int jmi = ks - 1;
		for (int i = k0; i < ks - 1; i++) {
			val = columnResponseExplanation(b, pvt, i, ks);
			if (val <= mi) {
				mi = val;
				jmi = i;
			}
		}
		return jmi;
	}

	/**
	 * Returns the squared ks-th response value if the j-th column becomes the
	 * ks-th after orthogonal transformation. A[][pvt[ks:j]] (or A[][pvt[j:ks]],
	 * if ks > j) and b[] are already QR-transformed on input and will remain
	 * unchanged on output.
	 * 
	 * More generally, it returns the inner product of the corresponding row
	 * vector of the response PaceMatrix. (To be implemented.)
	 * 
	 * @param b
	 *            PaceMatrix b
	 * @param pvt
	 *            pivoting vector
	 * @param j
	 *            the column A[pvt[j]][] is to be moved
	 * @param ks
	 *            the target column A[pvt[ks]][]
	 * @return the squared response value
	 */
	public double columnResponseExplanation(PaceMatrix b, IntVector pvt, int j,
			int ks) {
		/*
		 * Implementation:
		 * 
		 * If j == ks - 1, returns the squared ks-th response directly.
		 * 
		 * If j > ks -1, returns the ks-th response after
		 * Householder-transforming the j-th column and the response.
		 * 
		 * If j < ks - 1, returns the ks-th response after a sequence of Givens
		 * rotations starting from the j-th row.
		 */

		int k, l;
		double[] xxx = new double[n];
		int[] p = pvt.getArray();
		double val;

		if (j == ks - 1)
			val = b.A[j][0];
		else if (j > ks - 1) {
			int jm = Math.min(n - 1, j);
			DoubleVector u = getColumn(ks, jm, p[j]);
			DoubleVector v = b.getColumn(ks, jm, 0);
			val = v.innerProduct(u) / u.norm2();
		} else { // ks > j
			for (k = j + 1; k < ks; k++)
				// make a copy of A[j][]
				xxx[k] = A[j][p[k]];
			val = b.A[j][0];
			double[] cs;
			for (k = j + 1; k < ks; k++) {
				cs = g1(xxx[k], A[k][p[k]]);
				for (l = k + 1; l < ks; l++)
					xxx[l] = -cs[1] * xxx[l] + cs[0] * A[k][p[l]];
				val = -cs[1] * val + cs[0] * b.A[k][0];
			}
		}
		return val * val; // or inner product in later implementation???
	}

	/**
	 * QR transformation for a least squares problem<br/>
	 * A x = b<br/>
	 * implicitly both A and b are transformed. pvt.size() is the psuedo-rank of
	 * A.
	 * 
	 * @param b
	 *            PaceMatrix b
	 * @param pvt
	 *            pivoting vector
	 * @param k0
	 *            the first k0 columns of A (indexed by pvt) are pre-chosen.
	 *            (But subject to rank examination.)
	 * 
	 *            For example, the constant term may be reserved, in which case
	 *            k0 = 1.
	 **/
	public void lsqr(PaceMatrix b, IntVector pvt, int k0) {
		final double TINY = 1e-15;
		int[] p = pvt.getArray();
		int ks = 0; // psuedo-rank
		for (int j = 0; j < k0; j++)
			// k0 pre-chosen columns
			if (sum2(p[j], ks, m - 1, true) > TINY) { // large diagonal element
				steplsqr(b, pvt, ks, j, true);
				ks++;
			} else { // collinear column
				pvt.shiftToEnd(j);
				pvt.setSize(pvt.size() - 1);
				k0--;
				j--;
			}

		// initial QR transformation
		for (int j = k0; j < Math.min(pvt.size(), m); j++) {
			if (sum2(p[j], ks, m - 1, true) > TINY) {
				steplsqr(b, pvt, ks, j, true);
				ks++;
			} else { // collinear column
				pvt.shiftToEnd(j);
				pvt.setSize(pvt.size() - 1);
				j--;
			}
		}

		b.m = m = ks; // reset number of rows
		pvt.setSize(ks);
	}

	/**
	 * QR transformation for a least squares problem <br/>
	 * A x = b <br/>
	 * implicitly both A and b are transformed. pvt.size() is the psuedo-rank of
	 * A.
	 * 
	 * @param b
	 *            PaceMatrix b
	 * @param pvt
	 *            pivoting vector
	 * @param k0
	 *            the first k0 columns of A (indexed by pvt) are pre-chosen.
	 *            (But subject to rank examination.)
	 * 
	 *            For example, the constant term may be reserved, in which case
	 *            k0 = 1.
	 **/
	public void lsqrSelection(PaceMatrix b, IntVector pvt, int k0) {
		int numObs = m; // number of instances
		int numXs = pvt.size();

		lsqr(b, pvt, k0);

		if (numXs > 200 || numXs > numObs) { // too many columns.
			forward(b, pvt, k0);
		}
		backward(b, pvt, pvt.size(), k0);
	}

	/**
	 * Sets all diagonal elements to be positive (or nonnegative) without
	 * changing the least squares solution
	 * 
	 * @param Y
	 *            the response
	 * @param pvt
	 *            the pivoted column index
	 */
	public void positiveDiagonal(PaceMatrix Y, IntVector pvt) {

		int[] p = pvt.getArray();
		for (int i = 0; i < pvt.size(); i++) {
			if (A[i][p[i]] < 0.0) {
				for (int j = i; j < pvt.size(); j++)
					A[i][p[j]] = -A[i][p[j]];
				Y.A[i][0] = -Y.A[i][0];
			}
		}
	}

	/**
	 * Stepwise least squares QR-decomposition of the problem A x = b
	 * 
	 * @param b
	 *            PaceMatrix b
	 * @param pvt
	 *            pivoting vector
	 * @param ks
	 *            number of transformed columns
	 * @param j
	 *            pvt[j], the column to adjoin or delete
	 * @param adjoin
	 *            to adjoin if true; otherwise, to delete
	 */
	public void steplsqr(PaceMatrix b, IntVector pvt, int ks, int j,
			boolean adjoin) {
		final int kp = pvt.size(); // number of columns under consideration
		int[] p = pvt.getArray();

		if (adjoin) { // adjoining
			int pj = p[j];
			pvt.swap(ks, j);
			double dq[] = h1(pj, ks);
			int pk;
			for (int k = ks + 1; k < kp; k++) {
				pk = p[k];
				h2(pj, ks, dq[1], this, pk);
			}
			h2(pj, ks, dq[1], b, 0); // for matrix. ???
			A[ks][pj] = dq[0];
			for (int k = ks + 1; k < m; k++)
				A[k][pj] = 0;
		} else { // removing
			int pj = p[j];
			for (int i = j; i < ks - 1; i++)
				p[i] = p[i + 1];
			p[ks - 1] = pj;
			double[] cs;
			for (int i = j; i < ks - 1; i++) {
				cs = g1(A[i][p[i]], A[i + 1][p[i]]);
				for (int l = i; l < kp; l++)
					g2(cs, i, i + 1, p[l]);
				for (int l = 0; l < b.n; l++)
					b.g2(cs, i, i + 1, l);
			}
		}
	}

	/**
	 * Solves upper-triangular equation <br/>
	 * R x = b <br/>
	 * On output, the solution is stored in b
	 * 
	 * @param b
	 *            the response
	 * @param pvt
	 *            the pivoting vector
	 * @param kp
	 *            the number of the first columns involved
	 */
	public void rsolve(PaceMatrix b, IntVector pvt, int kp) {
		if (kp == 0)
			b.m = 0;
		int i, j, k;
		int[] p = pvt.getArray();
		double s;
		double[][] ba = b.getArray();
		for (k = 0; k < b.n; k++) {
			ba[kp - 1][k] /= A[kp - 1][p[kp - 1]];
			for (i = kp - 2; i >= 0; i--) {
				s = 0;
				for (j = i + 1; j < kp; j++)
					s += A[i][p[j]] * ba[j][k];
				ba[i][k] -= s;
				ba[i][k] /= A[i][p[i]];
			}
		}
		b.m = kp;
	}

	/**
	 * Returns a new matrix which binds two matrices together with rows.
	 * 
	 * @param b
	 *            the second matrix
	 * @return the combined matrix
	 */
	public PaceMatrix rbind(PaceMatrix b) {
		if (n != b.n)
			throw new IllegalArgumentException("unequal numbers of rows.");
		PaceMatrix c = new PaceMatrix(m + b.m, n);
		c.setMatrix(0, m - 1, 0, n - 1, this);
		c.setMatrix(m, m + b.m - 1, 0, n - 1, b);
		return c;
	}

	/**
	 * Returns a new matrix which binds two matrices with columns.
	 * 
	 * @param b
	 *            the second matrix
	 * @return the combined matrix
	 */
	public PaceMatrix cbind(PaceMatrix b) {
		if (m != b.m)
			throw new IllegalArgumentException("unequal numbers of rows: " + m
					+ " and " + b.m);
		PaceMatrix c = new PaceMatrix(m, n + b.n);
		c.setMatrix(0, m - 1, 0, n - 1, this);
		c.setMatrix(0, m - 1, n, n + b.n - 1, b);
		return c;
	}

	/**
	 * Solves the nonnegative linear squares problem. That is,
	 * <p>
	 * <center> min || A x - b||, subject to x >= 0. </center>
	 * <p>
	 * 
	 * For algorithm, refer to P161, Chapter 23 of C. L. Lawson and R. J. Hanson
	 * (1974). "Solving Least Squares Problems". Prentice-Hall.
	 * 
	 * @param b
	 *            the response
	 * @param pvt
	 *            vector storing pivoting column indices
	 * @return solution
	 */
	public DoubleVector nnls(PaceMatrix b, IntVector pvt) {
		int j, t, counter = 0, jm = -1, n = pvt.size();
		double ma, max, alpha, wj;
		int[] p = pvt.getArray();
		DoubleVector x = new DoubleVector(n);
		double[] xA = x.getArray();
		PaceMatrix z = new PaceMatrix(n, 1);
		PaceMatrix bt;

		// step 1
		int kp = 0; // #variables in the positive set P
		while (true) { // step 2
			if (++counter > 3 * n) // should never happen
				throw new RuntimeException("Does not converge");
			t = -1;
			max = 0.0;
			bt = new PaceMatrix(b.transpose());
			for (j = kp; j <= n - 1; j++) { // W = A' (b - A x)
				wj = bt.times(0, kp, m - 1, this, p[j]);
				if (wj > max) { // step 4
					max = wj;
					t = j;
				}
			}

			// step 3
			if (t == -1)
				break; // optimum achieved

			// step 5
			pvt.swap(kp, t); // move variable from set Z to set P
			kp++;
			xA[kp - 1] = 0;
			steplsqr(b, pvt, kp - 1, kp - 1, true);
			// step 6
			ma = 0;
			while (ma < 1.5) {
				for (j = 0; j <= kp - 1; j++)
					z.A[j][0] = b.A[j][0];
				rsolve(z, pvt, kp);
				ma = 2;
				jm = -1;
				for (j = 0; j <= kp - 1; j++) { // step 7, 8 and 9
					if (z.A[j][0] <= 0.0) { // alpha always between 0 and 1
						alpha = xA[j] / (xA[j] - z.A[j][0]);
						if (alpha < ma) {
							ma = alpha;
							jm = j;
						}
					}
				}
				if (ma > 1.5)
					for (j = 0; j <= kp - 1; j++)
						xA[j] = z.A[j][0]; // step 7
				else {
					for (j = kp - 1; j >= 0; j--) { // step 10
						// Modified to avoid round-off error (which seemingly
						// can cause infinite loop).
						if (j == jm) { // step 11
							xA[j] = 0.0;
							steplsqr(b, pvt, kp, j, false);
							kp--; // move variable from set P to set Z
						} else
							xA[j] += ma * (z.A[j][0] - xA[j]);
					}
				}
			}
		}
		x.setSize(kp);
		pvt.setSize(kp);
		return x;
	}

	/**
	 * Solves the nonnegative least squares problem with equality constraint.
	 * That is,
	 * <p>
	 * <center> min ||A x - b||, subject to x >= 0 and c x = d. </center>
	 * <p>
	 * 
	 * @param b
	 *            the response
	 * @param c
	 *            coeficients of equality constraints
	 * @param d
	 *            constants of equality constraints
	 * @param pvt
	 *            vector storing pivoting column indices
	 * @return the solution
	 */
	public DoubleVector nnlse(PaceMatrix b, PaceMatrix c, PaceMatrix d,
			IntVector pvt) {
		double eps = 1e-10 * Math.max(c.maxAbs(), d.maxAbs())
				/ Math.max(maxAbs(), b.maxAbs());

		PaceMatrix e = c.rbind(new PaceMatrix(times(eps)));
		PaceMatrix f = d.rbind(new PaceMatrix(b.times(eps)));

		return e.nnls(f, pvt);
	}

	/**
	 * Solves the nonnegative least squares problem with equality constraint.
	 * That is,
	 * <p>
	 * <center> min ||A x - b||, subject to x >= 0 and || x || = 1. </center>
	 * <p>
	 * 
	 * @param b
	 *            the response
	 * @param pvt
	 *            vector storing pivoting column indices
	 * @return the solution
	 */
	public DoubleVector nnlse1(PaceMatrix b, IntVector pvt) {
		PaceMatrix c = new PaceMatrix(1, n, 1);
		PaceMatrix d = new PaceMatrix(1, b.n, 1);

		return nnlse(b, c, d, pvt);
	}

	/**
	 * Generate matrix with standard-normally distributed random elements
	 * 
	 * @param m
	 *            Number of rows.
	 * @param n
	 *            Number of colums.
	 * @return An m-by-n matrix with random elements.
	 */
	public static Matrix randomNormal(int m, int n) {
		Random random = new Random();

		Matrix A = new Matrix(m, n);
		double[][] X = A.getArray();
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				X[i][j] = random.nextGaussian();
			}
		}
		return A;
	}

	/**
	 * Returns the revision string.
	 * 
	 * @return the revision
	 */
	public String getRevision() {
		return RevisionUtils.extract("$Revision: 1.6 $");
	}

	/**
	 * for testing only
	 * 
	 * @param args
	 *            the commandline arguments - ignored
	 */
	public static void main(String args[]) {
		System.out.println("================================================"
				+ "===========");
		System.out.println("To test the pace estimators of linear model\n"
				+ "coefficients.\n");

		double sd = 2; // standard deviation of the random error term
		int n = 200; // total number of observations
		double beta0 = 100; // intercept
		int k1 = 20; // number of coefficients of the first cluster
		double beta1 = 0; // coefficient value of the first cluster
		int k2 = 20; // number of coefficients of the second cluster
		double beta2 = 5; // coefficient value of the second cluster
		int k = 1 + k1 + k2;

		DoubleVector beta = new DoubleVector(1 + k1 + k2);
		beta.set(0, beta0);
		beta.set(1, k1, beta1);
		beta.set(k1 + 1, k1 + k2, beta2);

		System.out.println("The data set contains " + n + " observations plus "
				+ (k1 + k2)
				+ " variables.\n\nThe coefficients of the true model"
				+ " are:\n\n" + beta);

		System.out.println("\nThe standard deviation of the error term is "
				+ sd);

		System.out.println("==============================================="
				+ "============");

		PaceMatrix X = new PaceMatrix(n, k1 + k2 + 1);
		X.setMatrix(0, n - 1, 0, 0, 1);
		X.setMatrix(0, n - 1, 1, k1 + k2, random(n, k1 + k2));

		PaceMatrix Y = new PaceMatrix(X.times(new PaceMatrix(beta)).plusEquals(
				randomNormal(n, 1).times(sd)));

		IntVector pvt = (IntVector) IntVector.seq(0, k1 + k2);

		/*
		 * System.out.println(
		 * "The OLS estimate (by jama.Matrix.solve()) is:\n\n" + (new
		 * PaceMatrix(X.solve(Y))).getColumn(0) );
		 */

		X.lsqrSelection(Y, pvt, 1);
		X.positiveDiagonal(Y, pvt);

		PaceMatrix sol = (PaceMatrix) Y.clone();
		X.rsolve(sol, pvt, pvt.size());
		DoubleVector betaHat = sol.getColumn(0).unpivoting(pvt, k);
		System.out.println("\nThe OLS estimate (through lsqr()) is: \n\n"
				+ betaHat);

		System.out
				.println("\nQuadratic loss of the OLS estimate (||X b - X bHat||^2) = "
						+ (new PaceMatrix(X.times(new PaceMatrix(beta
								.minus(betaHat))))).getColumn(0).sum2());

		System.out.println("============================================="
				+ "==============");
		System.out.println("             *** Pace estimation *** \n");
		DoubleVector r = Y.getColumn(pvt.size(), n - 1, 0);
		double sde = Math.sqrt(r.sum2() / r.size());

		System.out.println("Estimated standard deviation = " + sde);

		DoubleVector aHat = Y.getColumn(0, pvt.size() - 1, 0).times(1. / sde);
		System.out.println("\naHat = \n" + aHat);

		System.out
				.println("\n========= Based on chi-square mixture ============");

		ChisqMixture d2 = new ChisqMixture();
		int method = MixtureDistribution.NNMMethod;
		DoubleVector AHat = aHat.square();
		d2.fit(AHat, method);
		System.out.println("\nEstimated mixing distribution is:\n" + d2);

		DoubleVector ATilde = d2.pace2(AHat);
		DoubleVector aTilde = ATilde.sqrt().times(aHat.sign());
		PaceMatrix YTilde = new PaceMatrix((new PaceMatrix(aTilde)).times(sde));
		X.rsolve(YTilde, pvt, pvt.size());
		DoubleVector betaTilde = YTilde.getColumn(0).unpivoting(pvt, k);
		System.out.println("\nThe pace2 estimate of coefficients = \n"
				+ betaTilde);
		System.out.println("Quadratic loss = "
				+ (new PaceMatrix(
						X.times(new PaceMatrix(beta.minus(betaTilde)))))
						.getColumn(0).sum2());

		ATilde = d2.pace4(AHat);
		aTilde = ATilde.sqrt().times(aHat.sign());
		YTilde = new PaceMatrix((new PaceMatrix(aTilde)).times(sde));
		X.rsolve(YTilde, pvt, pvt.size());
		betaTilde = YTilde.getColumn(0).unpivoting(pvt, k);
		System.out.println("\nThe pace4 estimate of coefficients = \n"
				+ betaTilde);
		System.out.println("Quadratic loss = "
				+ (new PaceMatrix(
						X.times(new PaceMatrix(beta.minus(betaTilde)))))
						.getColumn(0).sum2());

		ATilde = d2.pace6(AHat);
		aTilde = ATilde.sqrt().times(aHat.sign());
		YTilde = new PaceMatrix((new PaceMatrix(aTilde)).times(sde));
		X.rsolve(YTilde, pvt, pvt.size());
		betaTilde = YTilde.getColumn(0).unpivoting(pvt, k);
		System.out.println("\nThe pace6 estimate of coefficients = \n"
				+ betaTilde);
		System.out.println("Quadratic loss = "
				+ (new PaceMatrix(
						X.times(new PaceMatrix(beta.minus(betaTilde)))))
						.getColumn(0).sum2());

		System.out.println("\n========= Based on normal mixture ============");

		NormalMixture d = new NormalMixture();
		d.fit(aHat, method);
		System.out.println("\nEstimated mixing distribution is:\n" + d);

		aTilde = d.nestedEstimate(aHat);
		YTilde = new PaceMatrix((new PaceMatrix(aTilde)).times(sde));
		X.rsolve(YTilde, pvt, pvt.size());
		betaTilde = YTilde.getColumn(0).unpivoting(pvt, k);
		System.out.println("The nested estimate of coefficients = \n"
				+ betaTilde);
		System.out.println("Quadratic loss = "
				+ (new PaceMatrix(
						X.times(new PaceMatrix(beta.minus(betaTilde)))))
						.getColumn(0).sum2());

		aTilde = d.subsetEstimate(aHat);
		YTilde = new PaceMatrix((new PaceMatrix(aTilde)).times(sde));
		X.rsolve(YTilde, pvt, pvt.size());
		betaTilde = YTilde.getColumn(0).unpivoting(pvt, k);
		System.out.println("\nThe subset estimate of coefficients = \n"
				+ betaTilde);
		System.out.println("Quadratic loss = "
				+ (new PaceMatrix(
						X.times(new PaceMatrix(beta.minus(betaTilde)))))
						.getColumn(0).sum2());

		aTilde = d.empiricalBayesEstimate(aHat);
		YTilde = new PaceMatrix((new PaceMatrix(aTilde)).times(sde));
		X.rsolve(YTilde, pvt, pvt.size());
		betaTilde = YTilde.getColumn(0).unpivoting(pvt, k);
		System.out
				.println("\nThe empirical Bayes estimate of coefficients = \n"
						+ betaTilde);

		System.out.println("Quadratic loss = "
				+ (new PaceMatrix(
						X.times(new PaceMatrix(beta.minus(betaTilde)))))
						.getColumn(0).sum2());

	}
}
